home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / lin_alg.lha / lin_alg / matrix1.cc < prev    next >
C/C++ Source or Header  |  1993-08-08  |  13KB  |  551 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *            Linear Algebra Package
  6.  *
  7.  *        Basic linear algebra operations, level 1
  8.  *              Element-wise operations
  9.  *
  10.  ************************************************************************
  11.  */
  12.  
  13. #include "LinAlg.h"
  14. #include <math.h>
  15. #include <builtin.h>
  16.  
  17. #pragma implementation
  18.  
  19.  
  20.  
  21. /*
  22.  *------------------------------------------------------------------------
  23.  *            Constructors and destructors
  24.  */
  25.  
  26. void Matrix::allocate(
  27.     const short no_rows,        // No. of rows
  28.     const short no_cols,        // No. of cols
  29.     const short row_lwb = 1,    // Row index lower bound
  30.     const short col_lwb = 1        // Col index lower bound
  31. )
  32. {
  33.   valid_code = MATRIX_val_code;
  34.  
  35.   assure((nrows=no_rows) > 0, "No. of matrix cols has got to be positive");
  36.   assure((ncols=no_cols) > 0, "No. of matrix rows has got to be positive");
  37.  
  38.   Matrix::row_lwb = row_lwb;
  39.   Matrix::col_lwb = col_lwb;
  40.  
  41.   name = "";
  42.  
  43.   nelems = nrows * ncols;
  44.  
  45.   assert( (index  = calloc(ncols,sizeof(REAL *))) != 0 );
  46.   assert( (elements = calloc(nelems,sizeof(REAL))) != 0 );
  47.  
  48.   register int i;
  49.   for(i=0; i<ncols; i++)
  50.     index[i] = &elements[i*nrows];
  51. }
  52.  
  53.  
  54. Matrix::~Matrix(void)        // Dispose the Matrix struct
  55. {
  56.   is_valid();
  57.   delete index;
  58.   delete elements;
  59.   if( name[0] != '\0' )
  60.     delete name;
  61.   valid_code = 0;
  62. }
  63.  
  64.                 // Set a new Matrix name
  65. void Matrix::set_name(const char * new_name)
  66. {
  67.   if( name[0] != '\0' )
  68.     delete name;            // delete the old name (if any)
  69.   name = strcpy(calloc(strlen(new_name)+1,sizeof(char)),new_name);
  70. }
  71.  
  72.                 // Erase the old matrix and create a
  73.                 // new one according to new boundaries
  74.                 // with indexation starting at 1
  75. void Matrix::resize_to(const short nrows, const short ncols)
  76. {
  77.   is_valid();
  78.   if( nrows == Matrix::nrows && ncols == Matrix::ncols )
  79.     return;
  80.  
  81.   delete index;
  82.   delete elements;
  83.  
  84.   char * old_name = name;
  85.   allocate(nrows,ncols);
  86.   name = old_name;
  87. }
  88.                 // Erase the old matrix and create a
  89.                 // new one according to new boundaries
  90. void Matrix::resize_to(const short row_lwb, const short row_upb,
  91.                const short col_lwb, const short col_upb)
  92. {
  93.   is_valid();
  94.   const int new_nrows = row_upb-row_lwb+1;
  95.   const int new_ncols = col_upb-col_lwb+1;
  96.   Matrix::row_lwb = row_lwb;
  97.   Matrix::col_lwb = col_lwb;
  98.   if( new_nrows == Matrix::nrows && new_ncols == Matrix::ncols )
  99.     return;
  100.  
  101.   delete index;
  102.   delete elements;
  103.  
  104.   char * old_name = name;
  105.   allocate(new_nrows,new_ncols,row_lwb,col_lwb);
  106.   name = old_name;
  107. }
  108.  
  109.  
  110. /*
  111.  *------------------------------------------------------------------------
  112.  *             Making a matrix of a special kind    
  113.  */
  114.  
  115.                 // Make a unit matrix
  116.                 // (Matrix needn't be a square one)
  117. Matrix& Matrix::unit_matrix(void)
  118. {
  119.   is_valid();
  120.   register REAL *ep = elements;
  121.   register int i,j;
  122.  
  123.   for(i=0; i < nrows; i++)
  124.      for(j=0; j < ncols; j++)
  125.         *ep++ = ( i==j ? 1.0 : 0.0 );
  126.  
  127.   return *this;
  128. }
  129.  
  130.                 // Make a Hilbert matrix
  131.                 // Hilb[i,j] = 1/(i+j-1), i,j=1...max, OR
  132.                 // Hilb[i,j] = 1/(i+j+1), i,j=0...max-1
  133.                 // (Matrix needn't be a square one)
  134. Matrix& Matrix::hilbert_matrix(void)
  135. {
  136.   is_valid();
  137.   register REAL *ep = elements;
  138.   register int i,j;
  139.  
  140.   for(i=0; i < nrows; i++)
  141.      for(j=0; j < ncols; j++)
  142.         *ep++ = 1./(i+j+1);
  143.  
  144.   return *this;
  145. }
  146.  
  147. /*
  148.  *------------------------------------------------------------------------
  149.  *             Matrix-scalar arithmetics
  150.  */
  151.  
  152.                 // Assign a value to all the elements
  153. Matrix& Matrix::operator = (const double val)
  154. {
  155.   is_valid();
  156.   register REAL * ep;
  157.   for(ep=elements; ep < elements+nelems;)
  158.     *ep++ = val;
  159.  
  160.   return *this;
  161. }
  162.  
  163.                 // Add to all the elements
  164. Matrix& Matrix::operator += (const double val)
  165. {
  166.   is_valid();
  167.   register REAL * ep;
  168.   for(ep=elements; ep < elements+nelems;)
  169.     *ep++ += val;
  170.  
  171.   return *this;
  172. }
  173.  
  174.                 // Subtract a value from all the elements
  175. Matrix& Matrix::operator -= (const double val)
  176. {
  177.   is_valid();
  178.   register REAL * ep;
  179.   for(ep=elements; ep < elements+nelems;)
  180.     *ep++ -= val;
  181.  
  182.   return *this;
  183. }
  184.  
  185.                 // Multiply all the elements by a given value
  186. Matrix& Matrix::operator *= (const double val)
  187. {
  188.   is_valid();
  189.   register REAL * ep;
  190.   for(ep=elements; ep < elements+nelems;)
  191.     *ep++ *= val;
  192.  
  193.   return *this;
  194. }
  195.  
  196.                 // Return TRUE if all the elements are equal
  197.                 // to a given value
  198. int Matrix::operator == (const double val) const
  199. {
  200.   is_valid();
  201.   register REAL * ep;
  202.   for(ep=elements; ep < elements+nelems;)
  203.     if( *ep++ != val )
  204.       return 0;
  205.   return 1;
  206. }
  207.  
  208.                 // Return TRUE if all the elements are stricly
  209.                 // less than a given value
  210. int Matrix::operator < (const double val) const
  211. {
  212.   is_valid();
  213.   register REAL * ep;
  214.   for(ep=elements; ep < elements+nelems;)
  215.     if( !(*ep++ < val) )
  216.       return 0;
  217.   return 1;
  218. }
  219.  
  220.                 // Return TRUE if all the elements are stricly
  221.                 // greater than a given value
  222. int Matrix::operator > (const double val) const
  223. {
  224.   is_valid();
  225.   register REAL * ep;
  226.   for(ep=elements; ep < elements+nelems;)
  227.     if( !(*ep++ > val) )
  228.       return 0;
  229.   return 1;
  230. }
  231.  
  232. /*
  233.  *------------------------------------------------------------------------
  234.  *        Apply algebraic functions to all the matrix elements
  235.  */
  236.  
  237.                 // Take an absolute value of a matrix
  238. Matrix& Matrix::abs(void)
  239. {
  240.   is_valid();
  241.   register REAL * ep;
  242.   for(ep=elements; ep < elements+nelems; ep++)
  243.     *ep = ::abs(*ep);
  244.  
  245.   return *this;
  246. }
  247.  
  248.                 // Square each element
  249. Matrix& Matrix::sqr(void)
  250. {
  251.   is_valid();
  252.   register REAL * ep;
  253.   for(ep=elements; ep < elements+nelems; ep++)
  254.     *ep = *ep * *ep;
  255.  
  256.   return *this;
  257. }
  258.  
  259.                 // Take a square root of all the elements
  260. Matrix& Matrix::sqrt(void)
  261. {
  262.   is_valid();
  263.   register REAL * ep;
  264.   for(ep=elements; ep < elements+nelems; ep++)
  265.     if( *ep >= 0 )
  266.       *ep = ::sqrt(*ep);
  267.     else
  268.       info(),
  269.       _error("(%d,%d)-th element, %g, is negative. Can't take the square root",
  270.          (ep-elements) % nrows + row_lwb,
  271.          (ep-elements) / nrows + col_lwb, *ep );
  272.  
  273.   return *this;
  274. }
  275.  
  276.                 // Apply a user defined function
  277. Matrix& Matrix::user_function(USER_DEFINED_REAL_F uf)
  278. {
  279.   is_valid();
  280.   register REAL * ep;
  281.   for(ep=elements; ep < elements+nelems; ep++)
  282.     *ep = (*uf)(*ep);
  283.  
  284.   return *this;
  285. }
  286.  
  287. Matrix& Matrix::user_function(USER_DEFINED_DOUBLE_F uf)
  288. {
  289.   is_valid();
  290.   register REAL * ep;
  291.   for(ep=elements; ep < elements+nelems; ep++)
  292.     *ep = (*uf)(*ep);
  293.  
  294.   return *this;
  295. }
  296.  
  297. /*
  298.  *------------------------------------------------------------------------
  299.  *            Matrix-Matrix element-wise operations
  300.  */
  301.  
  302.                 // Check to see if two matrices are identical
  303. int operator == (const Matrix& im1, const Matrix& im2)
  304. {
  305.   are_compatible(im1,im2);
  306.   return (memcmp(im1.elements,im2.elements,im1.nelems*sizeof(REAL)) == 0);
  307. }
  308.  
  309.                 // Add the source to the target
  310. Matrix& operator += (Matrix& target, const Matrix& source)
  311. {
  312.   are_compatible(target,source);
  313.  
  314.   register REAL * sp = source.elements;
  315.   register REAL * tp = target.elements;
  316.   for(; tp < target.elements+target.nelems;)
  317.     *tp++ += *sp++;
  318.   
  319.   return target;
  320. }
  321.   
  322.                 // Subtract the source from the target
  323. Matrix& operator -= (Matrix& target, const Matrix& source)
  324. {
  325.   are_compatible(target,source);
  326.   register REAL * sp = source.elements;
  327.   register REAL * tp = target.elements;
  328.   for(; tp < target.elements+target.nelems;)
  329.     *tp++ -= *sp++;
  330.   
  331.   return target;
  332. }
  333.  
  334.                 // Modified addition
  335.                 //    Target += scalar*Source
  336. Matrix& add(Matrix& target, const double scalar,const Matrix& source)
  337. {
  338.   are_compatible(target,source);
  339.  
  340.   register REAL * sp = source.elements;
  341.   register REAL * tp = target.elements;
  342.   for(; tp < target.elements+target.nelems;)
  343.     *tp++ += scalar * *sp++;
  344.   
  345.   return target;
  346. }
  347.  
  348.                 // Multiply target by the source
  349.                 // element-by-element
  350. Matrix& element_mult(Matrix& target, const Matrix& source)
  351. {
  352.   are_compatible(target,source);
  353.   register REAL * sp = source.elements;
  354.   register REAL * tp = target.elements;
  355.   for(; tp < target.elements+target.nelems;)
  356.     *tp++ *= *sp++;
  357.   
  358.   return target;
  359. }
  360.  
  361.                 // Divide target by the source
  362.                 // element-by-element
  363. Matrix& element_div(Matrix& target, const Matrix& source)
  364. {
  365.   are_compatible(target,source);
  366.   register REAL * sp = source.elements;
  367.   register REAL * tp = target.elements;
  368.   for(; tp < target.elements+target.nelems;)
  369.     *tp++ /= *sp++;
  370.   
  371.   return target;
  372. }
  373.  
  374. /*
  375.  *------------------------------------------------------------------------
  376.  *            Compute matrix norms
  377.  */
  378.  
  379.                 // Row matrix norm
  380.                 // MAX{ SUM{ |M(i,j)|, over j}, over i}
  381.                 // The norm is induced by the infinity
  382.                 // vector norm
  383. double Matrix::row_norm(void) const
  384. {
  385.   is_valid();
  386.   register REAL * ep = elements;
  387.   register double norm = 0;
  388.  
  389.   while(ep < elements+nrows)        // Scan the matrix row-after-row
  390.   {
  391.     register int j;
  392.     register double sum = 0;
  393.     for(j=0; j<ncols; j++,ep+=nrows)    // Scan a row to compute the sum
  394.       sum += ::abs(*ep);
  395.     ep -= nelems - 1;            // Point ep to the beg of the next row
  396.     norm = ::max(norm,sum);
  397.   }
  398.   assert( ep == elements + nrows );
  399.  
  400.   return norm;
  401. }
  402.  
  403.                 // Column matrix norm
  404.                 // MAX{ SUM{ |M(i,j)|, over i}, over j}
  405.                 // The norm is induced by the 1.
  406.                 // vector norm
  407. double Matrix::col_norm(void) const
  408. {
  409.   is_valid();
  410.   register REAL * ep = elements;
  411.   register double norm = 0;
  412.  
  413.   while(ep < elements+nelems)        // Scan the matrix col-after-col
  414.   {                    // (i.e. in the natural order of elems)
  415.     register int i;
  416.     register double sum = 0;
  417.     for(i=0; i<nrows; i++)        // Scan a col to compute the sum
  418.       sum += ::abs(*ep++);
  419.     norm = ::max(norm,sum);
  420.   }
  421.   assert( ep == elements + nelems );
  422.  
  423.   return norm;
  424. }
  425.  
  426.  
  427.                 // Square of the Euclidian norm
  428.                 // SUM{ m(i,j)^2 }
  429. double Matrix::e2_norm(void) const
  430. {
  431.   is_valid();
  432.   register REAL * ep;
  433.   register double sum = 0;
  434.   for(ep=elements; ep < elements+nelems;)
  435.     sum += ::sqr(*ep++);
  436.  
  437.   return sum;
  438. }
  439.  
  440.                 // Square of the Euclidian norm of the
  441.                 // difference between two matrices
  442. double e2_norm(const Matrix& m1, const Matrix& m2)
  443. {
  444.   are_compatible(m1,m2);
  445.   register REAL * mp1 = m1.elements;
  446.   register REAL * mp2 = m2.elements;
  447.   register double sum = 0;
  448.   for(; mp1 < m1.elements+m1.nelems;)
  449.     sum += sqr(*mp1++ - *mp2++);
  450.  
  451.   return sum;
  452. }
  453.  
  454. /*
  455.  *------------------------------------------------------------------------
  456.  *             Some service operations
  457.  */
  458.  
  459. void Matrix::info(void) const    // Print some information about the matrix
  460. {
  461.   is_valid();
  462.   message("\nMatrix %d:%dx%d:%d '%s'",row_lwb,nrows+row_lwb-1,
  463.       col_lwb,ncols+col_lwb-1,name);
  464. }
  465.  
  466.                 // Print the Matrix as a table of elements
  467.                 // (zeros are printed as dots)
  468. void Matrix::print(const char * title) const
  469. {
  470.   is_valid();
  471.   message("\nMatrix %dx%d '%s' is as follows",nrows,ncols,title);
  472.  
  473.   const int cols_per_sheet = 6;
  474.   register int sheet_counter;
  475.   
  476.   for(sheet_counter=1; sheet_counter<=ncols; sheet_counter +=cols_per_sheet)
  477.   {
  478.     message("\n\n     |");
  479.     register int i,j;
  480.     for(j=sheet_counter; j<sheet_counter+cols_per_sheet && j<=ncols; j++)
  481.       message("   %6d  |",j+col_lwb-1);
  482.     message("\n%s\n",_Minuses);
  483.     for(i=1; i<=nrows; i++)
  484.     {
  485.       message("%4d |",i+row_lwb-1);
  486.       for(j=sheet_counter; j<sheet_counter+cols_per_sheet && j<=ncols; j++)
  487.     message("%11.4g  ",(*this)(i+row_lwb-1,j+col_lwb-1));
  488.       message("\n");
  489.     }
  490.   }
  491.   message("Done\n");
  492. }
  493.  
  494. #include <builtin.h>
  495. void compare(            // Compare the two Matrixs
  496.     const Matrix& matrix1,    // and print out the result of comparison
  497.     const Matrix& matrix2,
  498.     const char * title )
  499. {
  500.   register int i,j;
  501.  
  502.   are_compatible(matrix1,matrix2);
  503.  
  504.   message("\n\nComparison of two Matrices:\n\t%s",title);
  505.   matrix1.info();
  506.   matrix2.info();
  507.  
  508.   double norm1 = 0, norm2 = 0;        // Norm of the Matrixs
  509.   double ndiff = 0;            // Norm of the difference
  510.   int imax=0,jmax=0;            // For the elements that differ most
  511.   REAL difmax = -1;
  512.   register REAL *mp1 = matrix1.elements; // Matrix element pointers
  513.   register REAL *mp2 = matrix2.elements;
  514.  
  515.   for(j=0; j < matrix1.ncols; j++)    // Due to the column-wise arrangement,
  516.     for(i=0; i < matrix1.nrows; i++)    // the row index changes first
  517.     {
  518.       REAL mv1 = *mp1++;
  519.       REAL mv2 = *mp2++;
  520.       REAL diff = abs(mv1-mv2);
  521.  
  522.       if( diff > difmax )
  523.       {
  524.     difmax = diff;
  525.     imax = i;
  526.     jmax = j;
  527.       }
  528.       norm1 += abs(mv1);
  529.       norm2 += abs(mv2);
  530.       ndiff += abs(diff);
  531.     }
  532.  
  533.   imax += matrix1.row_lwb, jmax += matrix1.col_lwb;
  534.   message("\nMaximal discrepancy    \t\t%g",difmax);
  535.   message("\n   occured at the point\t\t(%d,%d)",imax,jmax);
  536.   const REAL mv1 = matrix1(imax,jmax);
  537.   const REAL mv2 = matrix2(imax,jmax);
  538.   message("\n Matrix 1 element is    \t\t%g",mv1);
  539.   message("\n Matrix 2 element is    \t\t%g",mv2);
  540.   message("\n Absolute error v2[i]-v1[i]\t\t%g",mv2-mv1);
  541.   message("\n Relative error\t\t\t\t%g\n",
  542.      (mv2-mv1)/max(abs(mv2+mv1)/2,1e-7) );
  543.  
  544.   message("\n||Matrix 1||   \t\t\t%g",norm1);
  545.   message("\n||Matrix 2||   \t\t\t%g",norm2);
  546.   message("\n||Matrix1-Matrix2||\t\t\t\t%g",ndiff);
  547.   message("\n||Matrix1-Matrix2||/sqrt(||Matrix1|| ||Matrix2||)\t%g\n\n",
  548.       ndiff/max( sqrt(norm1*norm2), 1e-7 )         );
  549.  
  550. }
  551.